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Abstract 


Large-eddy simulation results for laminar-to-turbulent transition in a spatially devel- 
oping boundary layer are presented. The disturbances are ingested into a laminar flow 
through an unsteady suction-and-blowing strip. The filtered, three-dimensional time- 
dependent Navier-Stokes equations are integrated numerically using spectral, high-order 
finite-difference, and three-stage low-storage Runge-Kutta methods. The buffer-domain 
technique is used for the outflow boundary condition. The localized dynamic model 
used to parameterize the subgrid-scale stresses begins to have a significant impact at 
the beginning of the nonlinear transition (or intermittency) region. The flow structures 
commonly found in experiments are also observed in the present simulation; the com- 
puted linear instability modes and secondary instability lambda-vortex structures are in 
agreement with the experiments, and the streak-like-structures and turbulent statistics 
compare with both the experiments and the theory. The physics captured in the present 
LES are consistent with the experiments and the full Navier-Stokes simulation (DNS), at 
a significant fraction of the DNS cost. A comparison of the results obtained with several 
SGS models shows that the localized model gives accurate results both in a statistical 
sense and in terms of predicting the dynamics of the energy-carrying eddies, without ad 
hoc adjustments. 


1. Introduction 

The problem of transition from laminar to turbulent flow in boundary layers is of great 
practical interest. Transition studies are motivated by a need to understand this physical 
process and to apply this knowledge to the prediction and control of transition. For example, 
the low skin-friction drag of laminar boundary layer flow compared with turbulent flow is 
very attractive to those who design high performance automobiles and aircrafts. On the 
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other hand, there are many cases where the high mixing and heat-transfer rates of turbulent 
boundary-layer are desirable, e.g., for combustion. 

Theoretical and experimental work to date in transition studies was reviewed in the recent 
paper by Kachanov (1994). It is well known today that natural transition due to the ampli- 
fication of small initial disturbances proceeds through a sequence of stages. The first stage is 
the two-dimensional (2-D) development of slowly growing (on viscous time-scales) Tollmien- 
Schlichting (TS) waves, which can be predicted accurately by the linear stability theory. The 
rapid growth (on convective time-scales) of infinitesimal three-dimensional disturbances on 
the base of finite-amplitude 2-D waves, which is explained by secondary instability theory, 
constitutes the second stage, and is followed by the appearance of small-scale motion and the 
final stages of transition. Herbert (1988) distinguished three types of secondary instability 
modes: fundamental, subharmonic and detuned. The fundamental and subharmonic instabil- 
ities lead respectively to aligned and staggered patterns of A-vortices observed in controlled 
experiments of transition. The detuned type leads to a combination pattern in which the 
A-vortices alternate between aligned and staggered patterns. If the initial disturbances are 
strong enough (for free-stream turbulence levels above 0.5%, Roach and Brierley 1992; Voke 
and Yang 1995), however, the transition will go into the nonlinear interaction stage directly, 
and the linear instability stage is bypassed; the transition is then called bypass transition. 
The present study focuses only on natural transition. 

It is now readily accepted that the initial stages, including receptivity and the linear and 
weakly nonlinear instability amplification, can be predicted quite accurately for many recep- 
tivity mechanisms or instabilities. What remains a significant challenge even with today’s 
supercomputers is the understanding and ability to compute and predict accurately the non- 
linear intermittency transition region. This ability is key for design because the maximum 
skin friction and temperature peaks occur in this region. Also, this region entails the presence 
of a great variety of scales to resolve. Hence, if a computational tool were available that could 
capture the relevant physics of the flow with only a modest cost penalty, it would provide 
the means to study the flow physics and validate proposed theories to model this complex 
region. 

Progress have also been made in the simulation capabilities of transition studies. Many 
of the previous numerical simulations were summarized in the review paper by Kleiser and 
Zang (1991). Although it is natural to simulate transition with a spatial approach, because 
transition in boundary-layer flows evolves in the streamwise direction, the majority of numer- 
ical simulations have involved computing the temporal growth of the instabilities. Spatial 
simulations of transition have been limited mainly by the extreme computational resolution 
requirements in the spatially-growing (streamwise) direction, and also by the problems associ- 
ated with assigning the proper inflow and outflow boundary conditions. It has been observed 
in experiments that the whole process of transition can require up to 25 wavelengths of the 
primary TS waves. Furthermore, the late stages of transition involve length scales between 
one and two orders of magnitude smaller than the TS scale. Thus, providing enough grid 
points to resolve this phenomena in such a long region is a daunting prospect. By contrast, 
in temporal simulations only one or two wavelengths need to be resolved, owing to the as- 
sumption of streamwise periodicity, and the assumption of parallel mean flow is generally 
made. 

Among the researchers that carried out spatial simulations of transition by direct simula- 
tions, Murdock (1986) simulated the K-type transition for the conditions of the Klebanoff et 
al. (1962) experiment, Fasel et al. (1990) simulated the early three-dimensional stages of both 
fundamental and sub harmonic transition. Recently, Bestek and coworkers (1994) calculated 
the 2D boundary layer transition under strong adverse pressure gradient, and Lundbladh 
et al. (1994) carried out the direct numerical simulation of bypass transition in both chan- 
nel and flat plate boundary layer flows. Joslin and Streett (1994) and Joslin (1995) have 
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simulated linear and nonlinear stationary crossflow vortex instability growth in swept wedge 
boundary-layer flows. Finally, Joslin (1995) used a fully 3D (non-periodic) DNS approach to 
simulate the evolution of instabilities in attachment-line flow. 

Large-eddy simulation (LES) is an alternative approach for the numerical solution of 
turbulent and transitional flow problems. Unlike the direct numerical simulation (DNS) in 
which all scales of motion are resolved, in LES only the dynamically important (large) scales 
are resolved, which results in significant decrease of the CPU time required for a simulation. 
The effect of the unresolved small scales on the large scales is modeled. Since the small 
scales tend to be more homogeneous and isotropic, the modelling can be simpler and more 
universal, compared with that in the Reynolds-averaged approach. 

While the application of LES to turbulent flows dates back to the seventies, only recently 
has this technique been used to the study of transitional flow. Piomelli et al. (1990), Piomelli 
and Zang (1991) and Germano et al. (1991) computed the transition in temporally developing 
boundary layer and plane channel flow. Their results indicate that, at the early stages of tran- 
sition, the eddy viscosity must be inactive to allow the correct growth of the perturbations. 
The dynamic model (Germano et al. 1991) achieves this result without the ad hoc corrections 
required by other models. Voke and Yang (1995) performed large-eddy simulation of bypass 
transition in a flat plate boundary layer, using a low-Reynolds-number correction for the 
Smagorinsky (1963) model; the properties of the simulated transition match those found ex- 
perimentally. Recently, Ducros et al. (1996) have performed an LES calculation of transition 
in a mildly compressible boundary layer using the Filtered Structure Function model. They 
obtain results that follow the expected trends in the early part of the transition region, but 
present no experimental or theoretical data to demonstrate quantitative agreement. In the 
late stages of transition and in turbulence, their resolution is too coarse, and the results, in 
the late breakdown stages and beyond, do not approach the expected turbulent laws. 

Large-eddy simulation, by definition, is a technique in which not all scales of motion are 
resolved. One question that may arise when applying it to transition problems regards the 
capability of LES to predict the development of shear layers and vortices whose scale is close 
the numerical filter. A subgrid-scale (SGS) model suitable for LES applications should not 
dissipate the energy of the low level perturbations during the initial stages of transition, but 
should reproduce the energy transfer to the unresolved scales during the nonlinear stages 
when such small, marginally resolved structures, are generated. 

Previous work on the application of LES to transitional flows concentrated on the pre- 
diction of statistical data (mean flow, Reynolds stresses, etc.) and only passing mention was 
devoted to the effect of the SGS model on the development of the vortical structures. This 
paper, on the other hand, will emphasize this issue, to illustrate how a localized eddy viscosity 
model responds to the local characteristics of the flow in such a way as to predict accurately 
the development of such structures, and compare this behavior with that of other models. 
In this work, we intend to simulate accurately the entire transition process, including the 
laminar, transitional and turbulent flow regimes. 

The numerical method will be introduced first, followed by the simulation description; 
the results will then be presented along with some comparisons with existing experimental 
and numerical data; the issue of the prediction of the vortical structures, and the effect of the 
eddy viscosity on the simulation will then be discussed. Finally, conclusions will be presented. 
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2. Problem formulation 


Mathematical formulation 

The present technique for simulating transitional flow relies on the decomposition of 
the instantaneous velocity fij(x, t) and pressure p(x,t) into base, f7 0 j(x) and P 0 (x), and 
disturbance, Uj(x, t) andp(x,t), components as 

fq(x,t) = Uoi{x.) T Uj (x, t) and p{x.,t) = P 0 (x) +p(x,i), (1) 


where x = («, y, z) are the streamwise, wall-normal, and spanwise coordinates and t is time. 
The base flow is the corresponding streamwise growing laminar flow field and is given by the 
Blasius similarity solutions. 

After substituting (1) into the Navier-Stokes equations and subtracting out the base-flow 
equations, the disturbance equations result and are given as 

duj duj duj dU oi _ dp 1 d 2 Uj 
dt J dx j 0J ' dxj 2 dxj dxj, Redxjdxj' 


and 


duj 

dx.: 


= 0 , 


( 3 ) 


where the Reynolds number Re = U 0 Q S* /u : U 0 c is the freestream velocity and S* is the 
boundary-layer displacement thickness at a reference streamwise location. The disturbance 
velocity boundary conditions are t/,;(x,t) = 0 at y = 0 and y — > oo. 

In LES, the large-scale (grid-resolved) components of the velocity and pressure are cal- 
culated and the effects of the small, unresolved scales are modeled. By applying the filtering 
operation 

7 ( x ) = / /(x’)G(x,x’)dx’ (4) 

JD 


(where G is the filter function and D is the entire domain) to (2) and (3), the governing 
equations for the large-scale velocity and pressure can be obtained: 


duj _ duj SfUi_ _ dU oi _ _&p_ _ drij_ 1 d 2 Uj 
dt ' <)./■ j ()x.j ' dx j dx, dxj Redxjdxj ’ 


and 



( 6 ) 


where r., ; is the subgrid-scale (SGS) stress tensor given by r, ; = uyuj — u,Uj. which must be 
modeled. 

Here, Tij is modeled by a localized version (Piomelli and Liu, 1995) of the dynamic SGS 
model (Germano et al. 1991). Local five point averaging in the streamwise and spanwise and 
three point averaging in the wall- normal direction is used in calculating the model coefficient. 
The total viscosity (molecular + eddy viscosity) is also forced to be non-negative to ensure 
numerical stability. 


Numerical method 

The filtered Navier-Stokes equations are solved using the fractional time step method 
(Chorin, 1968). Fourth-order finite difference and fourth-order compact difference schemes 
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Figure 1: Modulation functions for suction/blowing. 


are used in the streamwise direction (x) for the pressure and the momentum equations, 
respectively, Chebychev series are used in the wall- normal direction (y), and Fourier series in 
the spanwise direction ( z ). Implicit Crank-Nicolson time-advancement is used for the wall- 
normal diffusion terms and a three-stage Runge-Kutta scheme for the remaining terms. A 
code validation study was previously performed by Joslin et al. (1992, 1993). 

One of the major difficulties associated with the numerical simulation of spatially de- 
veloping boundary-layer transition is to define the outflow boundary conditions. Here the 
buffer domain technique proposed by Streett and Macaraeg (1989) is used, in which the gov- 
erning equations are gradually parabolized in a buffer region that is appended to the end 
of the computational domain, thus eliminating the necessity of applying outflow boundary 
conditions. 


Test-case parameters 

In the present simulation, the undisturbed laminar boundary layer is the base flow. Small 
disturbances are introduced into the flow through a suction/blowing strip at an upstream 
location. Since these disturbances are unstable, they are amplified as they propagate down- 
stream, and the flow goes through laminar, transition and turbulent stages consecutively. 

At the suction/blowing strip, a non-zero normal velocity v is applied. Two frequences are 
introduced into the flow by 

v(x,z,t) = Aif(x) sin (wt) + A 1/2 f(x) g(z) sin ( 7 ) 

where A\ and A\j 2 are the disturbance amplitudes, (f> is the phase shift between these two 
modes, and f{x), g(z) are the modulation functions shown in Fig. 1. The disturbance param- 
eters are chosen to match the experiment of Kachanov and Levchenko (1984) for controlled 
subharmonic breakdown and are given in Table 1. 

To simulate transition in a spatially developing boundary layer, one must use a computa- 
tional domain that is long enough in the streamwise direction so that laminar, transition and 
turbulent stages can all be observed, and the grid resolution must be fine enough to ensure 
accuracy. Two computational boxes were used for the present simulation (see Fig. 2). The 
first box covers the laminar region, the second one the transition and turbulent regions. In 
the first box, a fairly coarse grid can be used. In the second box, however, the grid resolu- 
tion has to be increased substantially to resolve the small scales present in the later stage 
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fundamental mode 

subharmonic mode 

frequency 

u) = 0.0916 

uj/2 

streamwise wavenumber 

a = 0.245 


streamwise wavelength 

A ;C = 25.6 


spanwise wavenumber 


/ 3 = 0.249 

spanwise wavelength 


A, = 25.3 

maximum amplitude 

A-i = 10 -3 

Ai = 10 -5 

2 

phase shift 

<f> = 0 

reference Re 

740 

inflow Re 

605 


Table 1: Parameters for the simulation of sub harmonic breakdown 


Box I 

Box II 

L x = 487, L y = 50, L z = 25 

L x = 359, L y = 100, L z = 25 

nx = 571, ny = 41, nz = 9 

nx = 561, ny = 65, nz = 33 

Ax = 0.86, Ay = 0.0085, A z = 3.2 

Ax = 0.64, Ay = 0.006, A z = 0.79 


Ax+ = 21, A y+ = 0.2, A^+ = 26 


Table 2: Dimensions of the computational boxes 


of transition. In both boxes, the grid points are evenly distributed in both streamwise (x) 
and spanwise ( z ) directions. In the wall-normal direction, the grid points are clustered in 
the near wall region. The box sizes, the grid sizes and the number of grid points used in the 
simulation are shown in Table 2, where Ay is the first wall normal grid size, i.e., the distance 
from the wall to the first grid point. This grid is finer by a factor of two in each direction 
than that used by Ductos et al. (1996). 

A buffer region with a streamwise length of is located at the end of each computational 
box for the numerical outflow boundary conditions. Numerical experiments indicate that this 
length is sufficient to allow the results in the useful region to be unaffected by the buffer region. 
The total lengths of the two computational domains are 19A_ r and 14Aj respectively. Both 
boxes have the same spanwise length A 2 . 

The suction/blowing strip is located 2A_ r downstream of the inlet and spans A ; „ and A 2 
in the streamwise and spanwise directions respectively. At the inflow of the first box, the 
disturbance velocity components ( u , v. w) are equal to zero. This is a valid assumption for 
the inflow boundary conditions because the suction/blowing strip is sufficiently far from the 
inflow. At the end of the first useful region (. x = 16A X ), the velocity data are saved every 
time step for two periods of the fundamental wave, i.e., one period of the sub harmonic wave. 
This information is used for the inflow boundary conditions for the second box. By using two 
periods of inflow data, no useful information is lost since waves with frequencies lower than 
that of the subharmonic mode are still not present in the flow up to this location. Numerical 
experiments indicate that, as long as the interface between the domains is located before 
the strongly nonlinear region, no significant disturbance is introduced by the interpolations 
required by this procedure. 
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Figure 2: Coordinate system and computational boxes. 


3. Simulation results 

In this section, the numerical results in the two useful regions of the computational boxes 
will be presented. In the first box, weakly nonlinear growth of the disturbances is observed. 
In the second box, the flow goes through laminar breakdown and the final transition to 
turbulence. Numerical results will first be shown, along with some comparisons with exper- 
imental and DNS data whenever possible, to validate the present results. The effects of the 
eddy viscosity on the resolved structures, as well as the intermittency model, will then be 
presented. 

Linear and weakly nonlinear stages 

In the first computational box, the fundamental wave initially grows before it reaches the 
second branch of the neutral curve, and the subharmonic modes grow much faster due to the 
secondary instability. Figure 3 shows the streamwise evolution of amplitude and phase of 
several harmonics for the streamwise velocity u , compared with the experiment by Kachanov 
and Levchenko (1984), where the Reynolds number is defined as Re x = xU^/v and x is the 
distance from the leading edge of the plate. Good overall agreement is obtained between the 
numerical simulation and the experiment. Towards the end of the computational domain, the 
simulation tends to give higher amplitudes for all the frequency modes shown here. In the 
experiment the fundamental mode (w) and its harmonics (2 uj and 3w) first saturate and then 
decrease, whereas in the simulation, their amplitude never decreases even though saturation 
can be seen. The same phenomena was observed in the direct numerical simulation of Fasel et 
al. (1990), who attributed it to insufficient resolution in the spanwise direction in their sim- 
ulation, because only the first and second harmonic components in z were included in the 
computation. However, in the present simulation the same trends are observed although the 
leading terms of nonlinear interactions with high harmonic components are included. For the 
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Figure 3: Streamwise evolution of amplitude (a) and phase (b) of selected harmonics at y/8 = 
0.26, z — 0 in the first computational box. Lines: LES; symbols: experiment (Kachanov and 

Levchenko, 1984). o : w; , •: co/2; , A : 2 u)\ , □ : 3w/2; , ■ : 

5w/2; , v 3lo. 
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same test case, Joslin et al. (1993) resolved the instabilities and raised the inconsistency issue 
between previous theoretical/computational comparisons and the Kachanov and Levchenko 
(1984) experiments. Imposing a nonzero pressure gradient on a variation in the recorded fun- 
damental frequency, they showed much better agreement in the computed and experimental 
fundamental modes; however, the subharmonic modes were in worse agreement. Their results 
suggest that there was insufficient information about the initial spectrum in the experiments 
to simulate adequately the same breakdown with computations. This could be due to the 
fact that in the simulation only two frequency modes are introduced into the flow, and non- 
linear wave interactions only occur between the fundamental mode, the subharmonic mode 
and other harmonics generated by nonlinear wave interactions, whereas in the experiment a 
broad spectrum of low-frequency modes, besides the subharmonic, was amplified to a signif- 
icant level from background noise, due to resonant interactions with the fundamental mode. 
Because of the existence and amplification of the low-frequency spectrum, wave interactions 
would then be much more complicated in the experiment, and transition in the experiment 
and in the LES may go through different routes; this is supported by the observation that 
the later stages of transition in the simulation feature high shear layer breakdown, whereas 
the experiment featured fast spectrum filling without spikes. It is also interesting to note 
that the fundamental wave and its harmonics tail up at the end of the useful region, which 
indicates the rapid joint growth of both fundamental and subharmonic waves at late stages 
of subharmonic breakdown after crossing their amplification curves, as was found in both 
experimental and theoretical studies (Maslennikova and Zelman, 1985; Crouch and Herbert, 
1993). 

The comparison of the streamwise phase evolution (Fig. 3b) shows excellent agreement 
between the simulation and the experimental data. All the modes have the same constant 
downstream propagation velocity throughout the whole box. This wave synchronism condi- 
tion, which should be satisfied in order to have resonant wave interactions between all these 
modes, was observed in the experiments of Corke and Mangano (1989) and Kachanov and 
Levchenko (1984). 

Wall-normal distributions of amplitude and phase of the fundamental and subharmonic 
modes for the velocity u are shown in Fig. 4. The comparison with the experiment is also 
good. 


Nonlinear transition stage and turbulent region 

In the second box, the disturbances grow further, all the frequency modes interact in a 
more complex manner, and the flow goes through laminar breakdown and final transition 
to turbulence. Fig. 5 shows the mean streamwise velocity in the wall- normal direction at 
five equally spaced streamwise locations, with the first one at the beginning and the last one 
at the end of the computational domain. The mean quantities shown in this section were 
obtained by averaging over both the spanwise direction and time, using 60 instantenous flow 
fields spanning 3 T in time (where T is the period of the fundamental wave) taken after a 
steady state was reached. As the flow evolves in the streamwise direction, the mean velocity 
changes from a laminar profile to a turbulent one. At the last location the velocity is in very 
good agreement with the DNS data of Spalart (1988), although the friction velocity for the 
LES is smaller than that for the DNS by 3%, due to the fact that the Reynolds number at 
this location is Reg* = 1372, higher than that of the DNS, Reg* = 1000. Figure 6 gives the 
rms of all three velocity components normalized by the friction velocity at the last streamwise 
location. Good agreement with DNS is again achieved. 

Figure 7 shows the streamwise variation of shape factor obtained from the large-eddy 
simulation. The shape factor in the turbulent region is slightly lower than that of the DNS, 
again due to Reynolds number difference. 
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Figure 4: Wall-normal distribution of amplitude (a) and phase (b) of the fundamental and 
subharmonic waves in the z = 0 plane, Re x = 3.6 x 10 5 . Lines: LES; symbols: experiment 
(Kachanov and Levchenko, 1984). , o : u>; , • co/2. 
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Figure 5: Mean streamwise velocity in wall units at various streamwise locations, +: DNS 
(Spalart, 1988); : U + = 2.44 In y + + 5.0; all other lines: LES. 



Figure 6: Wall-normal distribution of velocity rms at the last streamwise location, Re x = 
6.36 x 10 5 . Lines: LES, symbols: DNS. b : u rms ; ,*: v rms \ , o : w rms . 
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Figure 7: Streamwise variation of shape factor. : LES; : laminar flow; : DNS 


The local skin friction given in Figure 8 has the laminar value as long as the amplitudes 
of the disturbances are low, and begins to have a significant increase as the flow goes into 
laminar breakdown. The streamwise position where the skin friction begins to deviate from 
the laminar curve can be defined as the transition location, which in current calculation is 
Re x k, 4.6 x 10 5 . In the final stage of transition, the skin friction settles down to a level in 
agreement with the DNS (Spalart, 1988) of turbulent boundary layer at the corresponding 
Reynolds number Re§~ . The turbulent skin friction, denoted by the chain-dot-dot-dot line, 
is given by the equation (see White, 1991), 


0.455 

C/ “ ln 2 (0.06i?e iC ) ' 

with the virtual origin of the turbulent boundary layer moved from the leading edge towards 
the transition location. The turbulence skin friction thus obtained matches very well the 
result of LES. 


SGS model performance 

An important question for the simulation of transitional flows regards the capability of 
the LES to resolve the large vortical structures that arise during the nonlinear interaction 
stages of transition. In Fig. 9, the contours of the fluctuating streamwise velocity are shown 
in a horizontal plane located at y = 0.4 (y + = 13 in the turbulent region) away from the 
wall. The staggered A-vortices typical of the subharmonic breakdown can be first observed 
at x = 640 in the figure. They actually exist further upstream and by the time they are 
apparent in the figure, they have already become strong enough and gone into the later stage 
of development. The A-vortices are followed immediately by the laminar breakdown. Near 
the end of the computational box, the low speed streaks for the near-wall turbulent flow can 
be observed. 

Figure 10 shows a time-sequence of contours of the vertical shear du' /dy in the z = 0 
plane, where u! is the fluctuating velocity. A detached high shear layer first appears at 
x = 640 in Fig. 10(a), which is also the location where the skin friction and the shape factor 
start to depart from their laminar values; a kink develops on the shear layer at t = t 0 + 1.5T. 
Low pressure contour plots from the present simulation indicate a strong correspondence 
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Figure 8: Streamwise variation of local skin friction. : LES; : laminar flow 

0.664f?e -0 ' 5 ; : turbulent flow 0.455/ In 2 (0.06f?e_ r ) with the virtual origin at Re x = 

2.9 x fO 5 . 



Figure 9: Contours of fluctuating streamwise velocity v! in the x—z planes y = 0.38 (y + = f3) 
at t = t 0 + 1.5T, where t 0 = 756. The contour levels are between —0.15 and 0.15. 
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Figure 10: Time-sequence of vertical shear du 1 /dy in a x — y plane (2 = 0 ). (a) t = t 0 + T; 
(b) t = t 0 + 1.5T; (c) t = t 0 + 2 T; (d) t = t 0 + 2.5T, where i 0 = 756. The contour levels are 
between -0.45 and 0.45. 


between the kink in a shear layer and a vortex roll-up, as was also shown by Sandham and 
Kleiser (1992) by analyzing a temporal DNS database of the K-type breakdown. This small 
and energetic vortex can cause a spike in the velocity trace when it passes by the probe in 
an experiment. Multiple shear layer roll-ups occur in the subsequent development, and the 
flow finally reaches laminar breakdown and turbulence. In the turbulent region, the near- wall 
shear layers can be observed. It is worth pointing out that in the experiment of Kachanov 
and Levchenko (1984) no high-frequency spikes were observed, the transition featured a fast 
excitation of the broad spectrum of low-frequency fluctuations, including the subharmonic, 
followed by a filling of spectrum by an interaction of low-frequency fluctuations with the 
fundamental wave and its harmonics. 

In the present simulation, the A-vortices, the high shear layers and their roll-ups, and 
the streak structures are all captured, even though the grid resolution is marginal during the 
laminar breakdown, as indicated by the low level oscillations in the contours of Fig. 10. It 
can be concluded that even with much coarser grid than in DNS, LES using the localized 
dynamic model can predict transition accurately not only in terms of statistics, but even as 
far as the local behavior is concerned. 

One of the advantages of the dynamic model used in this calculation is its ability to adjust 
the model coefficient in time according to the flow conditions. By using the localized version 
of the dynamic model, the spatial variation of the flow can also be included. In the first box 
of the present simulation, the dynamic model gives essentially zero eddy viscosity because the 
transition is still in its early stages; this results in the correct prediction of the spatial growth 
of the disturbances. Figure 11 shows the wall-normal distribution of the mean eddy viscosity 
at the five evenly spaced streamwise locations shown in Fig. 5. At the first two streamwise 
locations, the eddy viscosity is essentially zero because the small scale disturbances are still 
not measurable in the flow; it becomes significant starting from the third location where 
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Figure 11: Mean eddy viscosity distributions at various streamwise locations. 


the laminar breakdown occurs; in the later stages of transition and in the turbulent stage it 
reaches the same magnitude as the molecular viscosity within the boundary layer. 

Figure 12 shows a time-sequence of contours of the SGS eddy viscosity in the z — 0 
plane. Comparing this figure with the corresponding Fig. 10 reveals again that the dynamic 
model only turns on the eddy viscosity in the transition region and afterwards, where small 
scales are present in the flow. Some negative values of eddy viscosity indicate backscatter 
(energy transfer from small scales to large scales) in the transition and the turbulent regions. 
Moreover, the figure suggests that the eddy viscosity is very small in the first spike stage, 
and becomes significant only during and after the multi-spike stage. This is clearly shown 
in Fig. 13 by the time-series of the disturbance velocity u and the SGS eddy viscosity at 
four locations in the flow field, whose positions were marked by the bullets in Figure 12a. 
In the laminar region where the disturbance is still small, the eddy viscosity given by the 
localized dynamic model is essentially zero (Fig. 13a). The first spike stage is indicated by a 
strong perturbation in the velocity trace at around to + 2.5 T, but the eddy viscosity remains 
small compared to the molecular one (Fig. 13b). In the multi-spike stage (Fig. 13c), the 
eddy viscosity has significant values at the time when multiple spikes appear in the velocity 
trace, because of the emergence of small scales; it, however, remains small before and after 
the spikes. Figure 13d shows the chaotic variation of the velocity and the eddy viscosity in 
the turbulent region. 

It has been shown that the localized SGS model used in the present study has the ability 
to adjust itself locally both in space and time according to the characteristics of the flow 
field. This is essential to the large-eddy simulation of transitional flows, since intermittency 
is a dominant phenomenon during transition. A SGS model that lacked this capability could 
result in over-damping of the disturbances and incorrect prediction of the transition process, 
as does the original Smagorinsky model. In addition, the eddy viscosity given by the localized 
model is zero in the first spike stage, and only multiple spikes can give rise to non-zero values 
of eddy viscosity. In a direct numerical simulation of transition, the grid resolution has to be 
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Figure 12: Time-sequence of the eddy viscosity in the (2 = 0) plane, (a) t = t 0 + T; (b) 
t = t 0 + 1.5T; (c) t = t 0 + 2T; (d) t = t 0 + 2.5T. The contour levels are between -1 and 3. 



increased before the multiple spike stage, since the effect of the small scale disturbances on 
the large ones then becomes significant and cannot be neglected. 

Also shown in Fig. 13 are the eddy viscosities obtained from simulations using the 
Smagorinsky model with the intermittency modification (Piomelli et al. 1990) 

v t = ^ ^ (Ac s ) \j2sijSij (9) 

where H is the shape factor, and the low-Reynolds-number model (Voke, 1995) 

Mr = v s ~ (^//?)[1 - exp {-(3v s /v)], (10) 

where v s is given by the original Smagorinsky model, c s is chosen as 0.1, and f3 = 4.5. 
According to Yoke (1996), only the fluctuating strain rate Sy, rather than the total strain, 
should be used to compute v s . The van Driest wall damping function is used for both models. 
The calculations start from the same flow field of LES at t = t 0 . It can be observed clearly 
that the eddy viscosity given by the low-Reynolds-number model is also insignificant in the 
laminar region, whereas in the first spike stage, the model gives much higher eddy viscosity; 
its variation closely resembles that of the vertical shear, reflecting the fact that du/dy is the 
dominant term in the strain-rate tensor; the spike is thus weakened. The higher viscosity 
of the low-Reynolds-number model also results in the smearing out of the multiple spikes in 
Fig. 13c. On the other hand, the Smagorinsky model with the intermittency modification 
gives results in agreement with the dynamic model well into the one-spike stage; it, however, 
gives smaller eddy viscosity in the multi-spike stage due to the intermittency attenuation. 
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Figure 13: Time history of the disturbance velocity u, the vertical shear du/dy and the eddy 
viscosity at y = 3.0, z = 0 and (a) x = 629; (b) x = 655; (c) x = 680; (d) x = 834. The 

eddy viscosities are given by: : localized dynamic model; : low-Reynolds-number 

Smagorinsky model; : Smagorinsky model with the intermittency modification. 
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Figure 14: Time-sequence of vertical shear du! jdy in a x — y plane (z = 0) and at t = t 0 + 2T. 
(a) Localized dynamic model; (b) low-Reynolds-number Smagorinsky model; (c) Smagorinsky 
model with the intermittency modification. The contour levels are between -0.45 and 0.45. 
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Figure 15: Streamwise variation of skin friction. : Localized dynamic model; : 

low-Reynolds number Smagorinsky model; : Smagorinsky model with the intermittency 

modification. 
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In the turbulent region, both models predict positive eddy viscosity and its peak values are 
smaller. 

Figure 14 compares the contours of du'/dy for different SGS models. The shear layer 
roll-up at x = 680 in the low-Reynolds-number model simulation is clearly not as strong as 
in other two cases, due to excessive damping provided by the eddy viscosity in the first spike 
stage; this also causes the transition to occur somewhat later and the increase in cj to be 
less steep, as shown in Fig. 15. The modified Smagorinsky model gives higher skin friction 
in the laminar breakdown region because of the attenuation of the SGS eddy viscosity by 
intermittency. The scattering of the data points in the turbulent region on Fig. 15 is due 
to insufficient averaging. The comparisons here between three models tends to indicate that 
the localized dynamic model does a better job, because the higher shear layer it resolves is 
sharper and the mean velocity profile it gives compares with DNS better than those of the 
other two. Furthermore, the localized dynamic model has no adjustable constant or ad hoc 
parameter. However, this issue need to be further investigated, a comparison with DNS is 
desired. 


4. Conclusions 

Large-eddy simulation of subharmonic type transition in a flat plate boundary layer was 
carried out. It was observed that the simulation predicts accurately the development of distur- 
bances in the early stages, the shear layers and vortical structures in the laminar breakdown 
stage, and the turbulent statistics in the turbulent stage. In addition, it has been shown that 
the skin friction for the simulation during laminar breakdown falls below the value of natu- 
ral transition. The localized dynamic SGS model has the desired capability to adjust itself 
locally both in space and time according to the characteristics of the flow field. It yields very 
small eddy viscosity in the laminar as well as the one spike stage, but gives significant values 
in the multi-spike and the turbulent stages, due to the emergence of energy in small scales 
of motion. The study also shows that the intermittency SGS model resembles the behavior 
of the localized dynamic model closely in the first-spike stage of transition, however, it gives 
smaller eddy viscosity in the multiple-spike stage; the low-Reynolds-number model predicts 
large values of eddy viscosity even in the first-spike stage. 
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